Dynamic Coupling - A New Tool for Three-Dimensional Flow and Transport Simulations

Simulating flow and transport in layered geologic media comprising aquifers and aquitards still represents a challenge to current finite-difference modeling techniques. In aquifers, the flow direction is usually horizontal and accurate simulation requires a sufficiently fine horizontal discretization. In aquitards, the flow direction is often vertical which requires grid refinement in the vertical direction. In such a situation, three-dimensional grids typically require tens of thousands of grid blocks. To address such problems, the ground-water simulator SWIFT 11 was enhanced to include a dynamic coupling option. When using this technique, the aquifers are discretized into normal finite-difference grids. However, the aquitards are discretized into columns of one-dimensional vertical subgrids, which are condensed from the main solution matrix. Dynamic coupling reduces the matrix of a full aquifer-aquitard system to the order and band width of the individual aquifer matrices, resulting in a significant reduction of memory requirements and execution time. Dynamic coupling was successfully applied in the context of developing and optimizing the aquifer remediation plan for a chemical-plant site in Tacoma, Washington. Past production of an arsenic-based herbicide has resulted in arsenic contamination of the groundwater beneath the plant. A remediation plan was developed to provide a longterm solution to the ground-water contamination. The enhanced version of SWIFT 11 was used to simulate three-dimensional ground-water flow and arsenic transport during and subsequent to the aquifer remediation. 92031° WMR WASTE 1 1 ~PDR INTRODUCTION Traditional finite-difference formulations transform the flow and transport equations into matrix equations. Most of the computer memory requirements and run-time demands are associated with storing and processing these matrices. Experience and theoretical studies have shown that memory requirements and computing time per time-step increase dramatically with the refinement of three-dimensional grids. Price and Coats (1973) have demonstrated that, using direct-solution algorithms with either standard ordering or more sophisticated ordering schemes for Gaussian elimination, the memory requirement can be calculated as follows: M = km n. (n. n2)2 (1) where nX, nY, nz = dimensions of the finite difference grid with n. 2 ny 2 nz km = ordering scheme specific constant (km = 1 for standard ordering). Standard ordering schemes number the nodes of a three-dimensional grid first along the shortest direction (z direction), then in the next shortest direction (y direction), and finally in the longest direction (x direction). The maximum band width of the matrix is then nX ny. Thus, the required memory increases with the square of the band width of the matrix. Price and Coats (1973) have further demonstrated that the computing time (work) increases with the cubed product of the two smaller grid dimensions: W = kw n, (nY n2) (2) where nx, nV. nz dimensions of the finite difference grid with n, ; n. 2 nz. kw ='ordering scheme specific constant (kw = 1 for standard ordering). Thus, the computing time increases with the cube of the band width of the matrix. Equations (1) and (2) imply that both M and W are most sensitive to increasing or decreasing the number of grid blocks along the shortest directions (i.e., y and z direction). In many applications, the z direction corresponds to the vertical direction of the model domain. Therefore, refining the model grid in the vertical direction often results in more than proportional increases of M and W. On the other hand, the largest gains are obtained from reducing the number of grid blocks along the vertical direction.